Spatial modelling etc

David L Miller & Mark V Bravington

International Whaling Commission Scientific Committee 2017

• Stratified Horvitz-Thompson is workhorse of many abundance estimates
• How is H-T going to fail?
• When do we need to use spatial models?
• What are “Bad surveys”?
• Report http://converged.yt/papers/iwc-2017-ht.pdf

• Today:
  1. what's wrong with H-T?
  2. Spatial models overview
  3. What can go wrong with spatial models?
  4. Testing designs in R
  5. Model checking for spatial models
• Tomorrow:
  1. What we missed, what is hard
  2. Try out tester on your data
  3. Other methods/software, future work
  4. Guidelines

• Try not to talk for more than an hour without a break
• We both have funny accents, yell if you don't understand!
• There is maths – don't worry

• This material usually takes 4 days+ to teach
• This will not prepare you to analyse spatial data
• BUT you can do this in St Andrews this summer!
• creem2.st-andrews.ac.uk

• Rescale the (flat) density and extrapolate

\[ \hat{N} = \frac{\text{study area}}{\text{covered area}}\sum_{i=1}^n \frac{s_i}{\hat{p}_i} \]

• \( s_i \) are group/cluster sizes
• \( \hat{p}_i \) is the detection probability (from distance sampling)

• Multiple sources of randomness in H-T equation:
  • \( \hat{p}_i \) - detectability
  • \( n \) - dealt with as \( n/L \), encounter rate
  • \( s \) - group size

• Probability of sampling every point in the study area is equal
• Is this true? Sometimes.
• If (and only if) the design is randomised

• Coverage probability
• H-T estimator assumes even coverage
• (or you can estimate)
• Otherwise not really valid

“A design is an algorithm for laying down samplers in the survey area”

“A realization (from that algorithm) is called a survey plan”

Len Thomas (Talk @CREEM 2004)

• Can't estimate w/ H-T w/o coverage
• “Fixed” “designs” violate assumptions
  • Some animals have \( \mathbb{P}(\text{included})=0 \)
• “Deteriorate” pooling robustness property
• What can we do?

• Encounter rate variance \( \approx n_j/l_j - n/L \)
• Within-transect variation can be bad
  • e.g., N-S transect, N-S density gradient

• If we suspect density change can stratify!
• Pre or post hoc (spatial and non-spatial)

• H-T is a spatial model (sort of)
• Violated an assumption if no randomness
  • Hard to assess how bad this is
• Fewster et al (2009) and Fewster (2011) give variance approaches

• Collect spatially referenced data
• Why not make spatially-explicit models?
• Go beyond stratified estimates
• Relate environmental covariates to counts

• Pure spatial, pure environmental, mixed?
• May have some prior knowledge
  • Biology/ecology
• What are drivers of distribution?
• Inferential aim
  • Abundance
  • Ecology

Hedley and Buckland (2004)

Miller et al. (2013)

Taking the previous example…

\[ \mathbb{E}\left(n_j\right) = \color{red}{A_j}\color{blue}{\hat{p}_j} \color{green}{\exp}\left[\color{grey}{ \beta_0 + \sum_k s_k(z_{kj})} \right] \]

\( n_j\sim \) some count distribution

• \( \color{red}{\text{area of segment}} \)
• \( \color{blue}{\text{probability of detection in segment}} \)
• \( \color{green}{\text{(inverse) link function}} \)
• \( \color{grey}{\text{model terms}} \)

• space, time, environmental (remotely sensed?) data

• 1-dimension: not much difference
• 2D more tricky
  • edge effects
  • tricky boundaries
  • more tomorrow
• Now going to do some maths…
  • (ignore at will)

• Fit needs to be penalised
• Something like:

\[ \int_\mathbb{R} \left( \frac{\partial^2 s(x)}{\partial x^2}\right)^2 \text{d}x\\ \]

• (Can always re-write this in the form \( \mathbf{\beta}^T \mathbf{S} \mathbf{\beta} \))
• Estimate the \( \beta_k \) terms but penalise objective
  • “closeness to data” + penalty (REML/ML)

• Generally:
  • penalties are improper prior precision matrices
  • (nullspace gives improper priors)
• Using shrinkage smoothers:
  
  • proper priors
  • empirical Bayes interpretation

library(dsm)
# environmental covariates
dsm_env_tw <- dsm(count~s(Depth) + s(NPP) + s(SST),
                  ddf.obj=df_hr,
                  segment.data=segs, observation.data=obs,
                  family=tw())
# space
dsm_xy_tw <- dsm(count~s(x, y), ddf.obj=df_hr,
                segment.data=segs, observation.data=obs,
                family=tw())

dsm is based on mgcv by Simon Wood

• Uncertainty from:
  • detection function parameters
  • spatial model
• Need to propagate uncertainty!
  • Methods in dsm
  • Bravington, Hedley & Miller (in prep)

• Build models in stages (detection function + GAM)
• Counts are functions of covariates
  • Pure spatial models
  • Environmental covariate models
  • Mix?!
• Fit/check using dsm
• Most of the theory is resolved, applications are hard

• Spatial modelling can give ubiased abundance ests
  • even with uneven coverage
  • limits to extrapolation
• V. even coverage => HT?
• “Evenness” subtle, detectability effect
  • e.g., weather bad in east

• Weather has a big effect on detectability
• Need to record during survey
• Disambiguate between distribution/detectability
• Potential confounding can be BAD

Thanks to Hiroto Murase and co for this data!

• Detection functions
  • Not enough observations
  • Uncertain species ID
  • Group size

Spatial models alone can't solve these issues

• Do you have enough observations?
• If they do look good (even coverage, etc)
• Is it worth re-analysing from H-T?
• Point estimates similar?
• Variance may well be different?

• Spatial models don't solve all problems
• Complex models can lead to complex issues
• Recording weather conditions is important
• You can always give up!

• Take a survey and simulate
• Is H-T robust?
• How do different spatial models compare?
  • Only thinking about total abundance & CV

• ltdesigntester github.com/dill/ltdesigntester
• (based on DSsim by Laura Marshall, CREEM)
• Setup simulations, test what can be done
• Most of the work needs to be done in GIS
  • Survey shapefiles, covariates etc
  • Import to R, runs models, shows output

• Spatial smoothers
  • thin plate spline, bs="tp" (Wood, 2003)
  • thin plate spline with shrinkage, bs="ts" (Marra et al., 2011)
  • Duchon spline, bs="ds", m=c(1, 0.5) (Miller et al., 2014)
  • tensor of thin plate spline (w/ and w/o rotated covariates)
• Stratified estimates
  • Horvitz-Thompson (w/ and w/o covariates)
  • stratified Horvitz-Thompson (w/ and w/o covariates)

• No model checking
• Dependent on “good” detection specs
• No group size model
• No \( g(0) \) or availability

library(ltdesigntester)
# setup a simulation
my_sim <- build_sim(design_path="path/to/shp",
                    dsurf=density_surface_matrix,
                    n_grid_x=dsurf_dim_x,
                    n_grid_y=dsurf_dim_y,
                    n_pop=true_N,
                    df=detection_function_specs,
                    region="path/to/shp")
# run it!
res <- do_sim(nsim=number_of_sims,
              scenario=my_sim,
              pred_dat=prediction_data_frame, ...)

• We can add covariates too (a wee bit unwieldy at the moment)
• Build multiple detection functions/sims in list()
• Covariates vary according to:
  • logit function E-W (can set pars, 2 state)
  • set values in segment data (already observed)

• We can test multiple (simple) scenarios
  • Assumption of simple gradients
  • Models likely won't work for difficult stuff if they don't work for simple things
• What will work/what won't
• Simple summary plots
• Better than the rest \( \neq \) good

• Count distribution
• Basis complexity
• Model (term) selection
• Sensitivity
• Observed vs. expected
• Cross-validation (replicability)

• Before: \( s(x) = \sum_{k=1}^K \beta_k b_k(x) \)
• How big should k be? “Big enough”
• Penalty takes care of the rest
• ?gam.check gives useful output
  • (also residual checks etc)

gam.check(dsm_env_tw)

Method: REML   Optimizer: outer newton
full convergence after 8 iterations.
Gradient range [-3.139726e-08,2.036272e-08]
(score 375.9503 & scale 4.316452).
Hessian positive definite, eigenvalue range [0.5725432,298.5906].
Model rank =  28 / 28 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

            k'   edf k-index p-value
s(Depth) 9.000 4.049   0.814    0.36
s(NPP)   9.000 2.846   0.779    0.04
s(SST)   9.000 4.916   0.771    0.04

“Everything is related to everything else, but near things are more related than distant things”

Tobler (1970)

• (approximate) \( p \) values (Marra & Wood, 2012)
  • path dependence issues
• shrinkage methods (Marra & Wood, 2011)
• ecological-level term selection
  • which biomass measure?
  • include spatial smooth or not?

• Diagnostic – compare observed vs. expected counts
  • Compare for different covariate/aggregations
• In next dsm, obs_exp() does this
• Going back to those rough POWER models…

> obs_exp(b, "beaufort")
               1        2       34
Observed 3.00000 10.00000 80.00000
Expected 6.97715 12.42649 83.03773

> obs_exp(b_nc, "beaufort")
                1        2       34
Observed 3.000000 10.00000 80.00000
Expected 8.478759 17.00705 73.23535

• How well does the model reproduce what we saw?
• Leave out one area, re-fit model, predict to new data
• Wenger & Olden (2012) have good spatial examples

• Not “cool” (statistically), but…
• Multi-stage models are handy!
• Understand and check each part
• Split your modelling efforts amongst people

• This methodology is general
  • Bears, birds, beer cans, Loch Ness monsters…
• Models are flexible!
  • Linear things, smooth things, random effect things (and more)
• If you know GLMs, you can get started with DSMs
  • Mature theoretical basis, still lots to do
• Active user community, active software development

Slides w/ references available at converged.yt

• For each \( b_k \) calculate the penalty
• Penalty is a function of \( \beta \)
  • \( \lambda \beta^\text{T}S\beta \)
• \( S \) calculated once
• smoothing parameter (\( \lambda \)) dictates influence

• We can set basis complexity or “size” (\( k \))
  • Maximum wigglyness
• Smooths have effective degrees of freedom (EDF)
• EDF < \( k \)
• Set \( k \) “large enough”

• “Fit to the histogram”
• Model:

\[ \mathbb{P} \left[ \text{animal detected } \vert \text{ animal at distance } y\right] = g(y;\boldsymbol{\theta}) \]

• Calculate the average probability of detection:

\[ \hat{p} = \frac{1}{w} \int_0^w g(y; \boldsymbol{\hat{\theta}}) \text{d}y \]

• Covariates that affect detectability (Marques et al, 2007)
• Perception bias (\( g(0)<1 \)) (Burt et al, 2014)
• Availability bias (Borchers et al, 2013)
• Detection function formulations (Miller and Thomas, 2015)
• Measurement error (Marques, 2004)  
  

Figure from Marques et al (2007)